#include "Hexa8NodeShpFunc.h"

using SG::DataStructure::Common::Real;
using SG::Algebra::Matrixd;
using SG::FEMSolver::Element::Hexa8NodeShpFunc;

Hexa8NodeShpFunc::Hexa8NodeShpFunc ()
:NSet (), DNDxiSet ()
{}

/** @brief 计算积分点处形函数及形函数对参数坐标导数
 * 
 * 缓存各个积分点处的计算结果，避免下一次积分点循环时重复计算（内存换取计算时间）
 */

void Hexa8NodeShpFunc::Compute (const SG::DataStructure::FEM::IntegPointArray& integPoints)
{
    const auto nGauss = integPoints.size ();
    NSet.resize (nGauss);
    DNDxiSet.resize (nGauss);
    const Real eighth{0.125};
    for (std::size_t index = 0; index < nGauss; ++index)
    {
        const Real& xi{integPoints[index].m_x};
        const Real& eta{integPoints[index].m_y};
        const Real& zeta{integPoints[index].m_z};

        NSet[index] = Matrixd (8, 1);
        NSet[index] (0, 0) = eighth * (1.0 - xi) * (1.0 - eta) * (1.0 - zeta);
        NSet[index] (1, 0) = eighth * (1.0 + xi) * (1.0 - eta) * (1.0 - zeta);
        NSet[index] (2, 0) = eighth * (1.0 + xi) * (1.0 + eta) * (1.0 - zeta);
        NSet[index] (3, 0) = eighth * (1.0 - xi) * (1.0 + eta) * (1.0 - zeta);
        NSet[index] (4, 0) = eighth * (1.0 - xi) * (1.0 - eta) * (1.0 + zeta);
        NSet[index] (5, 0) = eighth * (1.0 + xi) * (1.0 - eta) * (1.0 + zeta);
        NSet[index] (6, 0) = eighth * (1.0 + xi) * (1.0 + eta) * (1.0 + zeta);
        NSet[index] (7, 0) = eighth * (1.0 - xi) * (1.0 + eta) * (1.0 + zeta);

        DNDxiSet[index] = Matrixd (3, 8);
        DNDxiSet[index] (0, 0) = -0.125 * (1.0 - eta) * (1.0 - zeta);
        DNDxiSet[index] (1, 0) = -0.125 * (1.0 - xi) * (1.0 - zeta);
        DNDxiSet[index] (2, 0) = -0.125 * (1.0 - xi) * (1.0 - eta);

        DNDxiSet[index] (0, 1) =  0.125 * (1.0 - eta) * (1.0 - zeta);
        DNDxiSet[index] (1, 1) = -0.125 * (1.0 + xi) * (1.0 - zeta);
        DNDxiSet[index] (2, 1) = -0.125 * (1.0 + xi) * (1.0 - eta);

        DNDxiSet[index] (0, 2) =  0.125 * (1.0 + eta) * (1.0 - zeta);
        DNDxiSet[index] (1, 2) =  0.125 * (1.0 + xi) * (1.0 - zeta);
        DNDxiSet[index] (2, 2) = -0.125 * (1.0 + xi) * (1.0 + eta);

        DNDxiSet[index] (0, 3) = -0.125 * (1.0 + eta) * (1.0 - zeta);
        DNDxiSet[index] (1, 3) =  0.125 * (1.0 - xi) * (1.0 - zeta);
        DNDxiSet[index] (2, 3) = -0.125 * (1.0 - xi) * (1.0 + eta);

        DNDxiSet[index] (0, 4) = -0.125 * (1.0 - eta) * (1.0 + zeta);
        DNDxiSet[index] (1, 4) = -0.125 * (1.0 - xi) * (1.0 + zeta);
        DNDxiSet[index] (2, 4) =  0.125 * (1.0 - xi) * (1.0 - eta);

        DNDxiSet[index] (0, 5) =  0.125 * (1.0 - eta) * (1.0 + zeta);
        DNDxiSet[index] (1, 5) = -0.125 * (1.0 + xi) * (1.0 + zeta);
        DNDxiSet[index] (2, 5) =  0.125 * (1.0 + xi) * (1.0 - eta);

        DNDxiSet[index] (0, 6) =  0.125 * (1.0 + eta) * (1.0 + zeta);
        DNDxiSet[index] (1, 6) =  0.125 * (1.0 + xi) * (1.0 + zeta);
        DNDxiSet[index] (2, 6) =  0.125 * (1.0 + xi) * (1.0 + eta);

        DNDxiSet[index] (0, 7) = -0.125 * (1.0 + eta) * (1.0 + zeta);
        DNDxiSet[index] (1, 7) =  0.125 * (1.0 - xi) * (1.0 + zeta);
        DNDxiSet[index] (2, 7) =  0.125 * (1.0 - xi) * (1.0 + eta);
    }
}